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We investigate the dynamic heterogeneities of glassy particle systems in the theoretical schemes of 
bond breakage and four-point correlation functions. In the bond-breakage scheme, we introduce the 
structure factor Sh{q^ t) and the susceptibility Xh{t) to detect the spatial correlations of configuration 
changes. Here Xhif) attains a maximum at t = t^^^ as a function of time t, where the fraction of 
the particles with broken bonds 06 (t) is about 1/2. In the four-point scheme, treating the structure 
factor S^iq^t) and the susceptibility X4.{t), we detect superpositions of the heterogeneity of bond 
breakage and that of thermal low- frequency vibration modes. While the former grows slowly, the 
latter emerges quickly to exhibit complex space-time behavior. In two dimensions, the vibration 
modes extending over the system yield significant contributions to the four-point correlations, which 
depend on the system size logarithmically. A maximum of X4(^) is attained at t = t™^^, where these 
two contributions become of the same order. As a result, t^^^ is considerably shorter than t^^^. 

PACS numbers: 64.70. Q-,63.50.Lm,61.20.Lc,66.30.hh 



I. INTRODUCTION 

Recently, much attention has been paid to the dy- 
namics of glasses [1 . In particular, dynamic hetero- 
geneities exceeding the molecular size and emerging on 
long timescales |2J have been observed in a number of 
experiments [3]-[5] and molecular dynamics simulations 
in two dimensions (2D) and in three dimensions (3D) [5]- 
[T4l [T6H22] . In simulations, they can be detected if the 
spatial correlations of the particle configuration changes 
or the displacements between two separated times are cal- 
culated. In an early period, displacement heterogeneities 
were observed in applied strain in model amorphous al- 
loys [SHH]. Harrowell and coworkers visualized them in a 
one-component fluid [9 and a binary mixture [10^. Mu- 
ranaka and Hiwatari detected them on short [11] and 
long [T2 timescales in binary mixtures. Yamamoto and 
one of the present authors [13 -15 examined breakage of 
appropriately defined bonds and identified relatively ac- 
tive and inactive regions without and with applied shear 
flow. The bond-breakage events are produced by the con- 
figuration changes of the particle positions. The broken 
bonds accumulated in long time intervals are heteroge- 
neous such that their structure factor Sb{q,t) may be 
fitted to the Ornstein-Zernike form (ex 1/[1 + ^^^6(^)^]), 
where t is the interval width taken to be of the order of 
the structural relaxation time Tq,. The correlation length 
^5(t) grows with lowering the temperature T. Kob et 
al \16^ detected string-like motions of mobile particles 
as fundamental elements of structural relaxations, whose 
length distribution is widened with lowering T. 

Lacevic et a/. [23] presented a statistical theory of the 
dynamic heterogeneity in terms of the so-called four- 
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point dynamic correlation functions. They found that 
the four-point structure factor S4{q^t) can be fitted to 
the Ornstein-Zernike form and the susceptibility X4(^) 
exhibits a peak at a characteristic time ^4^^^ of order 
Tq,. The correlation length ^4 = ^4(^4^^^) thus obtained 
grows with lowering T. Subsequently, intensive efforts 
have been made to construct statistical theories and/or 
add further numerical results on the four-point correla- 
tions [iiiiHaQ]. 

However, there has been no systematic comparison 
between the bond-breakage scheme and the four-point 
scheme. The bond-breakage events occur as rare activa- 
tion processes, resulting in structural relaxations, in the 
absence of applied shear. In contrast, the physical pro- 
cesses giving rise to the four-point correlations have not 
yet been well understood. In this paper, we show that the 
four-point correlations originate twofold from the config- 
uration changes yielding the bond-breakage correlations 
and from the collective particle motions arising from the 
low- frequency transverse vibration modes [3TH4T] . The 
timescales of these two kinds of motions are dramati- 
cally different. In the latter, clusters of relatively mobile 
particles carry a large fraction of the vibrational energy 
and are distributed throughout the system [36]. The 
vibration modes have been studied to explain the low- 
temperature thermodynamic properties of glasses [1. 

In the low- frequency vibrational motions, the oscilla- 
tory particle displacements are highly heterogeneous so 
that the configuration changes should occur preferen- 
tially in more active regions with larger displacements, as 
pointed out by Schober et al [33 . This structural relax- 
ation mechanism was confirmed numerically in systems 
with particle numbers about 1000 [38l[39j and experimen- 
tally in quasi-2D colloidal glasses [42^. Thus, it explains 
the inseparable coupling between the structural disorder 
and the slow dynamics in glass. We mention some sim- 
ulations related to this coupling. Vollmayr-Lee et al [43] 
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found in 3D that mobile particles (in their definition) are 
surrounded by fewer neighbors than the others. Widmer- 
Cooper and Harrowell [44] detected a correlation between 
the short-time heterogeneity in a local Debye- Waller fac- 
tor and the long-time heterogeneity in 2D. Kawasaki et 
al. [22] claimed that medium-range crystalline order re- 
maining in glass controls ease of vitrification and nature 
of the glass transition. In polycrystal with small grains, 
the relation between the structure and the slow dynamics 
is more understandable, where the particles at the grain 
boundaries initiate configuration changes [45-47 . 

As a closely related effect, a very small applied strain 
produces strongly non-affine particle displacements in 
glass, indicating highly heterogeneous elastic moduli [6]- 
[Hl HSl |49]. Naturally, the particles in such elastically 
softer regions exhibit larger-amplitude displacements in 
the thermally excited vibration modes. Tanguy et a/. [35] 
showed that the classical elasticity theory holds only on 
spatial scales longer than a characteristic length (~30 
molecular sizes in their 2D model system). Moreover, in 
glass, irreversible plastic events are induced even by very 
small stains [M] [48] and plastic deformations are highly 
heterogeneous often leading to shear bands [5] |47l [50] . 
Under a fixed small strain at T = in 2D, Manning and 
Liu [51] numerically examined the relation between the 
low-frequency vibration modes and structural soft spots 
where configuration changes (particle rearrangements in 
their paper) are initiated. 

The organization of this paper is as follows. In Sec. II, 
our simulation method will be explained. In Sec. Ill, the 
bond-breakage scheme [T3[ [14] will be generalized. In 
Sec. IV, we will reexamine the four-point scheme [23] , 
where the collective particle motions arising from the 
vibartion modes will be identified. In Sec.V, the dy- 
namic heterogeneities detected by these two schemes will 
be compared. In Sec. VI, 3D results will be presented. 

II. NUMERICAL METHOD 

To illustrate consequences of the bond-breakage and 
four-point theories, we will show results of molecular dy- 
namics simulation of 50 : 50 binary mixtures composed 
of two species, 1 and 2, in 2D and 3D in amorphous 
states at low temperatures. We imposed the periodic 
boundary condition without applying shear fiow. The 
particle numbers of the two species are Ni = N2 = N/2. 
In 2D, TV win be mostly 4000 or 64000, but data for 
N = 16000 and 256000 wih also be given in Figs.3 and 8. 
In 3D, results for N = 10000 will be presented in Sec. VI. 
The two species have different diameters ai and (72 with 
o'2/o'i = 1.4 in 2D and (72/(11 = 1.2 in 3D. The particles 
interact via the soft-core potential, 

Va0ir)=e{^y^ -Ca^ (r<reut), (2.1) 

where a and P represent the particle species (a, /3 = 1,2), 
r is the particle distance, and e is the characteristic in- 



teraction energy. The interaction lengths are defined by 

^a/3 = (^a+^/3)/2. (2.2) 

The potential vanishes for r > rcut, where rcut = 4.5(7i 
in 2D and rcut = 3(7i in 3D. The constants Cai3 ensure 
the continuity of the potential at r = rcut- The masses 
of the two species satisfy 1712/ mi = ((72/(71)^. The aver- 
age number density is n = N/V = O.Sllcrf ^ in 2D and 
0.8(7]"^ in 3D, where V is the system volume. The system 
length L is 70.2(7i for N = 4000 and 281(7i for N = 64000 
in 2D, while L = 23.2(7i in 3D. Space and time will be 
measured in units of ai and 

^0 = (7i\/mi/e. (2.3) 

The temperature T will be measured in units of e/kB- 

We started from a liquid state at a high temperature, 
quenched the system to the final low temperature, and 
waited for a long time of order 10^. We imposed a ther- 
mostat in these steps. However, after this preparation 
of the initial states, we removed the artificial thermostat 
and integrated the Newton equations under the periodic 
boundary condition in the time range t > 0. This is 
needed to describe the effect of the vibration modes on 
long timescales. See the item (3) in the summary section 
for more discussions on the heat bath effect. Thus, the 
total particle number A/", the total volume and the 
total energy E are fixed in our simulation. 

III. BOND-BREAKAGE THEORY 

A. Background 

We regard two particles i and j with positions ri{t) 
and rj{t) to be bonded if [T3l[T4] 

rij{t) < Aiaai3, (3.1) 

where i e a and j G /3. Hereafter r^j(t) = \ri{t) — rj{t)\ 
is the distance between these particles at time t. At a 
later time t + At, this bond is treated to be broken if 

Tijit^At) > A2(7^^. (3.2) 

We assume that AiCJa^ is slightly larger than the peak 
distance of the pair correlation functions Qapij) and A2 is 
somewhat larger than A^. In this paper, we set Ai = 1.15 
and A2 = 1.5 in 2D and Ai = 1.3 and A2 = 1.7 in 3D. 

Let us consider the bonds at t = to and denote their 
total number as N}){to). A fraction of them will be bro- 
ken subsequently and the total number of the remaining 
bonds Ni){to + At) at t = to + At decays as 

Ni,{to + At)/iV5(to) = F5(At). (3.3) 

For large systems, the relaxation function ^^(At) may be 
treated to depend only on the time difference At (being 
nearly independent of the initial time to for large N). It 
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correlation function is expressed as 
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FIG. 1: Top: Bond relaxation function Fhit) in Eq.(3.3), self 
part of density time-correlation function Fs{q^ t) at q = 27r in 
Eq.(3.5), and fraction of non-S particles 1 — 06 (t) = 0s (t, 0) 
in Eq.(3.15) at T = 0.56 for N = 4000 in 2D. Relaxation 
times here are 8400 from Eq.(3.6), n 2.0 x 10^ ^ 35rc 
from Eq.(3.4). and Up = 4.18 x 10^ ^ Sr^ from Eq.(3.16). 
Bottom: fractions of B particles with k broken bonds 0s (t, k) 
in Eq.(3.14) for /c = 1, 2, and 3. 
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FIG. 2: Numerical results of (2nb/n)[l - Fb(t)], 0s(t, 1), and 
0s (t, 1) + 20s (t, 2) + 30s (t, 3) as functions of t in the early 
stage at T = 0.56 for N = 4000, which confirm Eq.(3.18). 



decreases with increasing At, so the bond-breakage time 
Tb may be defined by 

F,{n) = e-\ (3.4) 
On the other hand, the self part of the density time- 



= ^(J2exp[iq ■ Arj{to,to + t)] 



(3.5) 



where Arj{to,ti) = rj{ti) — rj{to) is the displacement 
vector of particle j and q is the wave vector. In Eq.(3.5), 
the average is taken over all the particles. In our simu- 
lation, the average over the initial time to and that over 
a number of runs were also taken. The structural relax- 
ation time Ta is usually defined at g = 27r by 



(3.6) 



Wave numbers will be measured in units of ^. 

In the upper panel of Fig.l, we display F})(t) and 
Fs{q,t) at (7 = 27r for T = 0.56 and TV = 4000, where 
To, = 8.4 X 10^ and n = 2.0 x 10^ = 35rc,. In the previous 
paper [14 , the relation = IOtc^ was found for binary 
mixtures with the soft-core potentential wih N = 10^ in 
3D. Both in 2D and 3D, Fb{t) may fairly be fitted to the 
stretched exponential form at low T as 



Fb{t)^exp[-{t/ny], 



(3.7) 



in the range t < ri^. At T = 0.56 we obtain c ~ 0.58 
in 2D. The exponent c increases with decreasing T. In 
contrast, Fs{q^t) exhibits a plateau /p(< 1) before the a 
relaxation at low T due to the thermal vibrational mo- 
tions (see the appendix). 



B. Bond-breakage correlations 

Here, we present a generalized formulation of the bond 
breakage to introduce broken-bond correlation functions. 
To this end, we define two overlap functions w^^l (r) and 



(2) 



(r) depending on the particle distance r as 



(3.8) 



with Ai and A2 being defined in Eqs.(3.1) and (3.2). The 
0{u) is the step function being equal to 1 for > and 
to for < 0. The fluctuating number density of the 
bonds may then be defined as 



(3.9) 



where we multiply 1/2 because a bond is supported by 
two particles in our definition. The statistical average of 
h}){r^t) is the average bond number density. 



If 1 



(3.10) 



Here, 715 ~ 3n in 2D at high densities. In fact, we numer- 
ically obtain = 2.28 = 2.81n for n = 0.811 in our 2D 
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system. Now we may introduce the broken bond number 
density in time interval [to^h] as 

V{r, to.h) = ^Yl ^^(^0' - ^^(^o)), (3.11) 

i 

where Bi{to,ti) is the broken bond number of particle i 
assuming a nonnegative integer quantity as 

Bdto,ti) = Ew^S(r.,(to))[l - wS(ri,(ii))]. (3.12) 

3 

This number tends to zero as ti to from Ai < A2 
and increases to 1, 2, • • • upon bond breakage. Hereafter, 
the particles with Si(to,ti) > 1 are called B particles, 
which are surrounded by different particle configurations 
at the initial and final times t = to and ti. On the other 
hand, those with Si(to, ti) = are called non-B particles, 
which have the same surrounding configurations at t = to 
and ti. The statistical average {V{r^to,ti)) depends on 
the time difference t = ti — to as 

Pb{t) = ^ J drV{r,toM) = ni,[l-F^{t)] (3.13) 

where is defined by Eq.(3.10) and Fi^it) by Eq.(3.6). 

Let the number of the particles with Si(to,ti) = k he 
Nsit, k) (A: = 0, 1, • • •) with t = ti - to. Then, 

Mt.k) = NB{t,k)/N (3.14) 

is the fraction of the B particles with k broken bonds for 
A: > 1, while (j)B{t, 0) = Nsit, 0)/N is the fraction of the 
non-B particles. The fraction of the total B particles is 
the following sum, 

Mt) = Y.Mt.k) = l- 0). (3.15) 

k>l 

We define the bond-preserving time r^p as 

1 - hiUp) = </>b(t6p,0) = e"^ (3.16) 

The particles have a broken bond on this time scale. 
Since each particle has several bonds (~ 6 in 2D), r^p is 
considerably shorter than the bond breakage time in 
Eq.(3.4). For t > ri^p^ the structural relaxation becomes 
appreciable. From Eq.(3.14) we obtain 

P6W = |E^'^^(*'^) (3.17) 

k 

From Eqs.(3.13) and (3.17) we find 

2'^[1 -F,{t)] = Y,kMt.k) (3.18) 

k 

Setting t = ti — to in steady states, we introduce the 
bond-breakage space-time correlation function, 

= ^(lZS^(^o,ti)Sfc(to,ti)(5(r-r,fc(to))y (3.19) 

^ ik ' 
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FIG. 3: Susceptibility Xh(f) in Eq.(3.21) vs t exhibiting a 
peak at t = for N = 16000. 

where Vi]^(to) = ri{to) — T/c(to). From Eq.(3.13) we have 
Gb{r^t) Pb{t)^ for large r. The structure factor of the 
broken bonds is given by 

Sb{q,t) = ^{\rq{to,ti)\^) 

= I dr[Gb{r,t)-pf,{tfy'i-'', (3.20) 

where q is the wave vector, q = \q\ is the wave num- 
ber, and Vq(to,ti) = (to, ^1) exp[zgr • 7-^(^0 )]/2 is 

the Fourier component of V{r,to,ti). 

As in the four-point scheme ^ , we introduce the sus- 
ceptibility Xb{t) to represent the overall degree of the 
bond-breakage correlations as follows: 

Xb{t) = ^(y25Bi{to,h)m{to,ti)), (3.21) 

^ ik ' 

in terms of the deviations 5Bi{to.>ti) = Bi{to^ti) — 
2pt{t)/n. Here (Bi) = T^^Bi/N = 2pb/n from Eqs.(3.11) 
and (3.13). In Sb{q^t) and Xb{t), the four particle posi- 
tions r^(to), rk{to), f'jiti)^ and r^(ti) are involved. In 
this sense, they are four-point correlation functions. 

Rabini et al. [52] introduced the number of particles 
that have left particle z's original neighbors at time t. It 
involves two times and was written as n°^*(0, t). It is sim- 
ilar to our Bi{to^to + 1) in Eq.(3.12). Abate and Durian 
[53] introduced a bond-breakage susceptibility Xb(^) sim- 
ilar to ours in Eq.(3.21) for a quasi- two-dimensional gran- 
ular system of air-fiuidized beads. 

C. Numerical results on bond breakage 

We further discuss consequences of our theory using 
numerical results in 2D. In Fig.l, we plot 1 — (l)b{t) = 
0b(^,O) in the upper panel and (psit^k) with /c = 1,2, 
and 3 in the lower panel, where r^p = 41800 = bra = 
0.14r5 from Eq.(3.16). At small (l)B{t^k) grow as 

(t)B{t,k) (X (3.22) 
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FIG. 4: (a) Structure factor 5'b(^, t™"^"") for bond breakage in 
Eq.(3.20) at various T for AT = 64000, (b) 5'6(^,t^^") at T = 
0.56 for various and (c) correlation lengths = ^bitb"^^) vs 
T for N = 64000. Sb{q,t) is approximately on a single curve 
independent of leading to weak system-size dependence of 
for K ^6 < L. 



where ai ~ 0.60, a2 ^ 1.0, and as ^ 1.3. Though we 
cannot derive these exponents theoretically, they should 
arise from correlated occurrence of bond breakage events. 

In Fig. 2, we confirm the validity of Eq.(3.18) from 
simulation at T = 0.56, where rib/n = 3.17. We find 
that (2n5/n)[l — Fi){t)] nearly coincides with 1) for 

t < 200 and with 1) + 20b(^,2) + 3^b(^,3) for 

t < 4000 within a few % differences. The first relation 
for k = 1 in Eq.(3.22) is consistent with Eq.(3.18) since 
ai = c. In fact, for t <^ri), Eqs.(3.7) and (3.18) yield 



(3.23) 



In Fig. 3, we plot Xb{t) in Eq.(3.21) as a function of t for 
N = 16000, which exhibits a peak at t = t^^^. Here, we 
have ^5(t) ~ 0.5 at t = t^^^, which can be seen in Fig.l 
for N = 4000 and T = 0.56. Here, however, Xb{t) begins 
to increase for t > t^^^, because the B particles with 
Bi = 2,3, ••• becomes appreaciable at very long times 
(see Fig.l). In (a) of Fig. 4, we show Sb{q,t^^^) vs q for 
various T with N = 64000. In its calculation we took the 
average over the initial time to in a wide range of [0, 10^] 
for T > 0.64 and that of [0, 2 x 10^] for T = 0.56, which 
was needed because of the slow bond breakage. We may 
fairly fit Sb{q^t) to the Ornstein-Zernike form [l3l[T4], 



S,{q,t) = x',{t)/[l^Q'm' 



(3.24) 



where Xbi^) ~ l™g^o >^b{Q^ t) is the long wavelength limit 
of Sb{q^t) and ^5 = ^b{t) is the correlation length rep- 



resenting the spatial scale of the correlated configura- 
tion changes. Furthermore, in (b), we show Sb{q^t^^^) 
vs q at T = 0.56 for various which demonstrates 
weak system-size dependence of the bond-breakage cor- 
relations. For N = 256000, however, the averaging over 
the initial time to is still insufficient because of very 
large t^^^ ~ 4 x 10^ as compared to the simulation time 
(~ 10^). As a result, the corresponding S^iq^t^^^) ex- 
hibit noticeable fluctuations at small q. In (c), the corre- 
lation length ^b vs T is plotted at t = t^^^ for TV = 64000, 
which increases with lowering T. Note that ^5 = ^5(^5^^^) 
is nearly independent of the system size from (b) as long 
as K ^5 < L. 

In the original papers [I3l [M], the broken-bond struc- 
ture factor was defined differently, so it is written as 
Sj^iq^t) here. It was calculated for the Fourier com- 
ponent of the following broken bond number density. 



^Yo(r,to,ti) = ^E^S(^^^-(^o))[l-^S(n,(ti))] 



xS{r - Rij(to)), 



(3.25) 



where t = ti — to was set equal to 0.05r5 or O.lrfo. Here, 
the midpoint position Rij(to) = |(^i(^o) -^'^j{to)) of the 
two particles i and j is used instead of the position ri{to) 
in S{r — ri{to)) in Eq.(3.11). Note that the particles 
pairs with a common broken bond are included in Sb{q^ t). 
As a result, the inter-particle peak at g = 27r appears 
in Sb{q,t) (which is not shown in Fig. 4), while it does 
not apppear in Sj^{q^t). However, there is no essential 
difference between these two definitions for q < 2. In 
addition, in the previous work [14 , the dynamic scaling 



relation of the form ' 
in 2D and z = 2 in 3D. 



was obtained, where z = 4 



IV. FOUR-POINT THEORY 

Lacevic et al. [23] introduced the four-point correlation 
function to analyze the dynamic heterogeneity in glassy 
systems. In their numerical analysis of a 3D binary mix- 
ture in the NVE ensemble, they used the Lennard- Jones 
potential, where the particle size ratio was (72 /ai = 1.2 
and the particle numbers were Ni = N2 = 4000. We crit- 
ically review their theory comparing it with our theory 
of bond breakage using some numerical analysis. 



A. Overlap and nonoverlap with initial regions 

For a time interval [to?^i] (^ = ^1 — ^0 > 0), we intro- 
duce a fluctuating density variable, 

Q(r, to, ti) = ^ ^z(to, h)S{r - r,{to)). (4.1) 

i 

For each i we define a nonnegative integer, 

Tiito^ti) = Y.w{\ri{to) - r,(ti)|), (4.2) 
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FIG. 5: Time evolution of Ari(to,to + t) for five parti- 
cles among 4000 particles. They exhibit thermal fluctua- 
tions with magnitudes on the order of the overlap length 
0.3. Two particles escape from the initial circle at t ^ 8000 
(blue A) and ~ 10^ (red A). They have a broken bond with 
Bi{to,to + t) = 1 and Ti{to,to + t) = in later times. For 
the other three particles, the bonds are preserved, but fre- 
quently fluctuate between 1 and in the time range displayed. 
Inset: Time evolution of the time-smoothed distance Afi{t) 
defined in Eqs.(4.8) and (4.10) in the range 4000 < t < 8000. 



using the following overlap function [23] 
w{r) — 6{A4ai — r). 



(4.3) 



The overlap length A/^ai is common for the two parti- 
cle species for simplicity. In Eq.(4.2) the particle posi- 
tions r^(to) and rj{ti) are those at different times. Thus 
J~'i{to^h) is the number of overlapping particles in the 
initial circle (or sphere in 3D) \r — ri{to)\ < A^ai in two 
configurations separated by time t = ti — to. We may 
call Q(r,to,^i) the two-point overlap density. 

In the original analysis [23 , the overlap function was 
written as w{r) = 6{aa2 — r) with a(j2 = 1.2acri, where 
the parameter a was set equal to 0.3 maximizing the 
four-point susceptibility X4^(f) (see the discussion around 
Eq.(4.17) below). In numerical analysis in this paper, we 
set = 0.3 in Eq.(4.3). These selected values are con- 
siderably shorter than the particle radii, but somewhat 
exceed the square root of the plateau value of the mean 
square displacement [23 . Thus, as ti ^ to, the distinct 
terms with j ^ i vanish in the summation of Eq.(4.2), 
leading to J^i{to,ti) 1 and Q{r,to,ti) p(r,to), 
where p{r^t) = ^^S{r — ri{t)) is the usual fluctuating 
number density. For ti > to, J^i{to^ti) is either or 
1. We did not detect the particles with > 2 in our 
simulation. 

In the definition of Q(7",to,ti), the terms in Eqs.(4.1) 
and (4.2) may be divided into the self part with i = j 
and the distinct part with i j. The self part of Q reads 

Qs{r, to, h) = ^ w{Au{to,h))S{r - r,(to)), (4.4) 




FIG. 6: (a) Trajectory of the particle escaping from a cage 
(red line in Fig.5). Data points at t = 8000 + 20A: {k = 
0, 1, • • • , 180) are written. Its position is (x, y) = (23.36, 67.00) 
at t = 10140 (point A) and {x,y) = (24,44,67.16) at t = 
10160 (point B). Between these times, the particle escapes 
from the circle \r — ri{to)\ < 0.3 (red line), (b) Displacements 
between this short time interval [10140, 10160]. Several par- 
ticles undergo a string-like motion, where the particle (A) in 
the upper panel has two broken bonds and the others have 
one broken bond. 



where Ari(to,ti) is the displacement length of particle z, 



Ar,(to,ti) = |r,(ti)-n(to)|. 



(4.5) 



Setting a = 0.3, Lacevic et a/. [23] found that the self part 
gives rise to the dominant contributions in the four-point 
correlations. Also in our simulation, the self part dom- 
inates over the distinct part. We consider the particles 
with Ari > 0.3 and = 1, for which another particle j 
has moved within the initial circle of particle i giving rise 
to the distinct contribution. We compare their number 
with the number of the particles with Ar^ > 0.3. For ex- 
ample, in a 2D simulation run for T = 0.56 and N = 4000 
(which yielded Fig. 13 and the left panels of 14), these 
numbers are 256 and 2195, respectively, at t = 10200 
(see Table II). In a 3D simulation run for T = 0.56 and 
TV = 4000 (which yielded Figs. 15 and 16), they are 366 
and 1452, respectively, at t = 10^ (see Table IV). More 
than 90% of these distinct particles (with Ar^ > 0.3 and 
J^i = 1) are B particles having broken bonds. 
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(a) At^200 (b)At^400 




FIG. 7: Early-stage snapshots of particles with J^^ = at t = 
200 (left) and 400 (right), which are classified into those with 
i3i = (red) and those with Bi = 1 (blue). Here T = 0.56 and 
N — 4000. Numbers of the former (red) and the latter (blue) 
are (766, 75) in (a) and (535, 81) in (b). The former arise from 
the collective vibrational modes emerging for t > 10. 



B. Background vibrational fluctuations 

We first examine the fluctuations of the particle posi- 
tions due to the thermally excited low-frequency vibra- 
tion modes at low T before the onset of the structural re- 
laxation. They exhibit significant heterogeneity, as found 
by Muranaka and Hiwatari in a very short time interval 
of width 5 in a 2D soft-core system with N = 10^ [11 . 

For t <^T})p^ we may neglect the configuration changes 
from the discussion below Eq.(3.16) and define the dis- 
palcement vector Ui(t) = ri(t) — fi^ where fi is the 
time-averaged position in a time interval with width 
much shorter than r^p. The equal-time variance {\u\^) = 
^i{\ui\^) / N is a half of the plateau value Mp of the mean 
square displacement at low T (see Eq.(A5) below). From 
the appendix, {\u\^) increases with increasing the system 
length L logarithmically in 2D as [54] 

= Mp/2^Co + Ciln(i/<Ji), (4.6) 

where Co and Ci are functions of T and are independent 
of L. In terms of the shear modulus ji and and the bulk 
modulus the coefficient Ci is expressed as 

In our 2D system, we obtained /i = 18 and K = 67.5 
at T = 0.56 from the inital linear growth of stress- strain 
relation and the density-pressure relation (not shown in 
this paper) [55^. Then we find Ci ^ 0.0060 at T = 0.56. 
If the total particle number is increased from 4000 to 
64000, Eq.(4.6) yields the incremental increase Ci ln4 = 
0.0083 in In fair agreement with this estimate, we 

numerically calculated to be 0.02266 for N = 4000 

and 0.03305 for N = 64000, where f ^ was equated with 
the time-average of ri{t) in a time interval of width 200. 

With Eq.(4.6), we need to examine how the particles 
remain within or go outside their initial circle. In Fig. 5, 
we display Ari{to^ti) in Eq.(4.5) for five particles in time 



range t = ti - to < 18000 with N = 4000. These parti- 
cles are separated from one another with distances longer 
than 10. We recognize that these displacements undergo 
rapid thermal fluctuations with magnitudes nearly equal 
to the overlap length 0.3. In the early stage {t < ri,p)^ 
most of J^i fluctuate between 1 and 0. In this example, 
two of them escape from their initial circle, where one 
has a broken bond at t ~ 8000 and the other has two 
broken bonds dit t ^ 10000 (see Fig. 6). Each jump it- 
self occurs on a short time of order 10. In (a) of Fig. 6 
gives the trajectory of the particle escaping from a cage 
at t ~ 10^ in Fig. 5. In (b) of Fig. 6, this escape take place 
as a string-like motion involving several particles as in 
3D [T6l[33]. See (a) of Fig. 13 below for other examples. 
In (b) of Fig. 6, we can see that most of the particles in- 
volved have only one broken bond, which is particularly 
the case for isolated configuration changes. Therefore, 
these motions may also be treated as small slips in 2D. 

Removing the rapid temporal fluctuations, we calcu- 
lated the smoothed displacement lengths, 

Afi{t) = \fi{t^to)-ri{to)l (4.8) 
for the time-smoothed positions, 

1 rtsm 

riit) = T- / dt'nit + t'). (4.9) 

^sm Jo 

In the inset of Fig. 5, the smoothing time tsm is 500. Even 
on this timescale, the particles move considerably, even 
across their initial circle. Note that this tsm is much 
longer than the traversal time of the transverse sounds 
across the system ta = L/c± ~ 17 for A/" = 4000 [55] , 
These complex fluctuations should originate from su- 
perposition of weakly coupled, low-frequency vibration 
modes |3iH4T]. 

In Fig. 7, we show the particles with = at very 
early times t = 200 and 400 for T = 0.56 and TV = 4000, 
where their initial circles contain no particle. The frac- 
tion of the particles with = 0, written as 04 (t), is soon 
about 0.2 for t > 10. That is, a considerable amount of 
the non-B particles with = already appear from 
very early times. However, some of the B particles de- 
picted at t = 200 are changed to those with = 1 at 
t = 400 (see the sentences at the end of Sec.IVA). 

In Fig. 8, we show the particles with = at t = 
10"^ ~ Ta for T = 0.56 in a much larger system of 
N = 64000. Some heterogeneities have sizes of order 
50. Among the displayed particles, the B and non-B 
particles amount to 22034 and 18122, respectively. The 
patterns of the latter are more extended than those of 
the former. Moreover, their timescales are distinctly sep- 
arated (see Fig. 14 below). These indicate the presence of 
thermally excited large-scale vibration modes. 

C. Four-point correlations 

The average of Q(r,to,ti) in Eq.(4.1) is written as 
q4{t) = (Q(r,to,ti)) 



FIG. 8: Late-stage snapshot of particles with J^^ = in a 
large system of = 64000 and L = 281, where t = 10^ and 
T = 0.56. Depicted particles are classified into those with 
Bi — (red) and those with Bi > 1 (blue). Numbers of 
the former (red) and the latter (blue) are (22034, 18122), re- 
spectively. The former arise from the low-frequency vibration 
modes. Subsequent time-evolution in the upper box region 
will be given in the right of Fig. 14. 
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FIG. 10: (a) Four-point structure factor S^iq^t^^^) vs q in 
Eq.(4.13) for various T, where N = 64000. (b) S^iq^tT") 
vs q for various iV at T = 0.56. (c) ^A{tT") for N = 64000. 
Marked system-size dependence appears at small q due to the 
low-frequency vibration modes. 




FIG. 9: Susceptibility x^ifi) in Eq.(4.14) vs t for iV = 4000. 



w{\ri{h)-rj{h)\), (4.10) 



which is a function oi t = ti — t^. In our simulation, 
= 1 or 0, so (74 (t) is related to the fraction 04 (t) of the 
particles with = as 



(4.11) 



Here, ^4(0) = n and ^4(00) = vqih? ^ where vq is the area 
or volume of the overlap region {vq = ^{A^ai)'^ in 2D 
and vo = 47r(A4cri)^/3 in 3D). 

As in Eq.(3.19) for Ghir^t)^ the four-point space-time 
correlation function is given by 

G4(r,t) = {Q{r^r'MM)Q{r'MM)) 



UY.J^^{to.h)J^k{to.tl)S{r - r,k{to))). (4.12) 



V 



where rik{to) = Ti(to) — Tfe(^o)- The four-point structure 
factor is defined by 

S4{q,t) = ^{\Qq{to,ti)\^) 

= J dr[G,{r,t)-qi{tfy'i-'', (4.13) 

where Qq(to,*i) = J^j{to,h) exp[iq ■ rj{to)] is the 
Fourier component of Q{r,to,ti). We define the four- 
point susceptibihty Xiit) by 



Xiit) 



V 



Y,SJ'i{to,h)5J'k{to,ti)). (4.14) 



ik 



in terms of the deviation STi{to, ti) = Ti{to^ ti) — q4{t)/n. 
In these correlation functions, the four particle positions 
^^(^o), f'kito)^ f'jiti)^ and r£{ti) are involved. However, 
as discussed around Eq.(4.4), the self parts with i = j 
and k = £ dominate over the distinct parts with i ^ j 
andk^i [23 . 

Berthier et al.[24l[25j showed that the four-point sus- 
ceptibility X4(t) depends on the ensemble {NVE or 
NVT) and the dynamics (Newtonian or Brownian). Our 
X4(^4^^^) in the NVE ensemble is roughly 60% of the 
long wavelength limit of the four-point structure factor 
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xlitT'''') = ^^^q^o S^iq^tf^"^), which is consistent with 
the previous calculations [24H28] . However, in our cal- 
culation, the long wavelength limit of the bond-breakage 
structure factor Xbi'^b^^^) — l™g^o 'S'5((7, tJJ^^^) cannot be 
determined reliably because of the very long t^^^ at low 
T (see Fig. 4) and our maximum bond-breakage suscep- 
tibility X6(^6^^^) apparently exceeds Xh^'^T^^) 
a few ten %. Thus, we cannot draw a definite conclusion 
on the relation between Xbl^b^^"") ^^^d Xh^^h"^^)- 

In (a) of Fig. 9, we give X4(^) Eq.(4.14) as a func- 
tion of t for various T, which is calculated in the NVE 
ensemble with N = 4000. It is maximized at t = ^4^^^. 
Here, due to the transverse sound propagation, a smaller 
acoustic peak emerges with lowering Tatt = ta/2~8.6 
[5F, whose existence has not been reported in the previ- 
ous papers. It becomes more evident at lower T, where 
the acoustic damping is weaker. We also calculated X4(^) 
for other N. For = 64000, the acoustic peak was at 
t ~ 34 and its height even exceeded the first peak height 
for low T. For TV = 1000, there was no acoustic peak. 
See the item (3) in the summary for more discussions. 

Lacevic et al. [23] determined the overlap length a(j2 to 
maximize the peak height of the four-point susceptibility 
X4(^4^^^) as a function of the parameter a at T = 0.59 
in 3D. They then obtained a = 0.3 and used it also at 
other low T. Following their method, we also maximized 
X4(^4^^^) as a function of the overlap length to obtain 
= 0.3 in Eq.(4.3) for T = 0.64 and TV = 4000. Then 



t] 



10"^ at T = 0.56 both for TV 



4000 and 64000. 
In (a) of Fig. 10, we plot S4{q^t^^^) vs q for various T 
with TV = 64000. In its calculation, we took the average 
over the initial time to in the wide range [0, 10^]. As in 
the case of Sb{q^t) in Eq.(3.20), we may fit S4,{q^t) to the 
Ornstein-Zernike form as E3| 



(4.15) 



where X4(^) = lii^ig^o 5*4 ^) is the long wavelength limit 
of S/^i^q^t) and £,^{1) is the four-point correlation length. 
Furthermore, in (b), we show S4,{q^f^^^) at T = 0.56 
for various TV up to TV = 256000 to demonstrate its sig- 
nificant system-size dependence at small q. In (c), we 
display ^4(^4^^^) vs T for TV = 64000 as an example, 
which increases with lowering T. We recognize that the 
ratio i^it^f^^) / (.hi^h"^^) exceeds unity and increases with 
increasing TV. For example, it is about 3 for T = 0.56 
and TV = 64000. Thus, on our 2D simulation, the effect 
of the low-frequency vibration modes on the four-point 
correlations becomes stronger for larger TV. 



COMPARISON OF THE TWO 
THEORETICAL SCHEMES 




FIG. 11: Fraction 06 (t) and susceptibility Xh{i) for bond 
breakage and those 04 (t) and X4(t) for four-point correlations 
for T = 0.64 and TV = 4000. Here, Xb{t) and X4(t) are maxi- 
mized for (j)h{t) ^ 1/2 and 04 (t) ^ 1/2, respectively. 
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FIG. 12: Relaxation times as functions of T for the soft-core 
potential for TV = 4000 in 2D. From above, they are Th in 
Eq.(3.4), Up in Eq.(3.16), t^^^ from the maximum of Xh{t), 
^max £j,Qj^ maximum of x^if)^ and in Eq.(3.6). 



and 04 (t) is that with J^i = 0. We can see that 05 ~ 1/2 
at t = t^^"" and 04 1/2 at t = tf^"". If 05 (t) (or 04 (t)) 
is close to or 1, Si){q^ t) (or S\){q^ t)) becomes very small. 

So far we have introduced the bond-breakage time r5 in 
Eq.(3.4), the relaxation time Tq, from Fs{q^t) in Eq.(3.6), 
the bond-preserving time r^p in Eq.(3.16), the maximiza- 
tion time t^^^ of Xb{t) in (a) of Fig. 3, and the maximiza- 
tion time ^4^^^ of X4:{t) in Fig. 9. In Fig. 12, these times 
are in the following order, 



A. Timescales 

In Fig. 11, we display 05 (t), 04 (t), Xb{t), and X4.{t) at 
T = 0.64 for TV = 4000 in the two theoreticaal schemes. 
Here, 05 (t) is the fraction of the particles with Bi > 



' Up < n, 



(5.1) 



with 35 < n/Ta < 10^ in the range 0.56 < T < 0.96 
for TV = 4000. In our 2D simulation, and tf^"" 
exhibit strong system-size dependence due to the low- 
frequency vibration modes. In fact, we numerically ob- 
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(a) An 



{b)8. 




•••• . 



5«o 




FIG. 13: Snapshots at t = ti - to = 10^ for T = 0.56 and N = 4000 in 2D. System length is L = 70.2. (a) Displacements 
Ari(to,ti) with \ Ari\ > 0.3, (b) B particles, and (c) those with = classified into B particles (blue) and non-S particles 
(red). Non-S particles with = are produced by collective motions, while most B particles participate in string-like motions 
and satisfy J^j = 0. 




FIG. 14: Snapshots at t = 10000 + 200A: {k = 1,2) at T = 0.56 for N = 4000 (left) and for N = 64000 (right). Top: 
Particles with = classified into non-S particles (red) and B particles (blue), where the former significantly change but 
the latter little change in this short time interval. The change of the former is larger for the larger N in the right. See the 
preceding snapshots in (c) of Fig. 13 for N — 4000 and in the black box of Fig. 8 for N — 64000. Bottom: Displacements 
An = Ari{to + 1', to + + 200) multiplied by 5 with | Ar^| > 0.3 in short time intervals of 200, where t' = 10000 or 10200. 



tained = 8400 for N = 4000 and = 2140 for 
N = 64000 (see the appendix). It is worth noting that 
significant system-size dependence of the plateau behav- 
ior of Fs{q,t) (and Tq, from Eq.(3.6)) has been reported 
in 2D and 3D simulations [271 IMl |57] . Furthermore, Kar- 



makar et al. [27| examined system size-dependence of ^4 
and S^iq.t) up to N = 351232 in 3D. 

In addition, we comment on the stress-time-correlation 
function, which considerably decreases in the early stage 
due to the thermal motions as well as Fs{q^t). As a 
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TABLE I: Particle numbers with (a) Bi > 0, (b) Bi > and 
Ti = 0, (c) Bi = = 0, and (d) An > 0.3 in Fig. 14 at 
t = 10200 and 10400 for N = 4000 and 64000 in 2D. Those 
of the particles common in these two shots are also given. 
Subscript i is omitted from An, Bi, and J^i. 



t{N) 


B>0 


B>0,J- = 


B = J^ = 


Ar > 0.3 


10200(4000) 


1629 


1104 


1091 


2451 


10400(4000) 


1667 


1180 


1221 


2643 


common(4000) 


1576 


925 


823 


1923 


10200(64000) 


26188 


18336 


20884 


43544 


10400(64000) 


26413 


18678 


21785 


44693 


common(64000) 


25192 


14149 


14349 


31790 



TABLE IL Numbers of non-B and B particles for three cat- 
egories in Fig. 14 at t = 10200 for N = 4000 in 2D. 





Ar > 0.3 

= 


Ar > 0.3 
1 


Ar < 0.3 
1 


total 


B = 


1091 


7 


1273 


2371 


B>0 


1104 


249 


276 


1629 


total 


2195 


256 


1549 


4000 



result, its relaxation time is of order [30^, while the 
nonlinear flow behavior is characterized by [M]. 



B. Time-evolution on long and short timescales 

In (a) of Fig. 13, the arrows represent relatively large 
displacements Avi = ri{ti) — ri{to) with lAr^l > 0.3 
[TTJ[T8]. We can see both large- amplitude string-like mo- 
tions and smaller-amplitude collective motions. In (b), 
all the B particles are displayed. In (c), those with 
J^i = are divided into B and non-B particles, where 
the former exhibit patterns closely resembling those of 
the B particles in (b). The number of the B particles 
with J^i = is about 70% of that of the total B particles 
in (b). 

From (a) and (c) of Fig. 13, the non-B particles with 
J^i = mostly arise from the collective motions on large 
scales, as has been the case at short times in Fig. 7. Their 
selection is very sensitive to the overlap length, aa2 in 
the original work [23 and A^ai in this paper, while the 
B particles (even with J^i = 0) are relatively insensitive 
to it. In the four-point theory |23], the overlap length 
was chosen to maximize xT^^ — XA.i^'T^^)^ discussed 
above Eq.(4.15). Roughly speaking, their method is to 
maximize the contribution from the thermal collective 
motions of non-B particles to XA^if) (see Fig. 11). 

We next examine time-evolution at two consecutive 
times t = 10200 and 10400 for N = 4000 (left) and 64000 
(right) at T = 0.56. The four upper panels of Fig. 14 dis- 
play the particles with J-i(to, to -\- 1) = grouped into B 
and non-B particles. These snapshots are subsequent to 
that in (c) of Fig. 13 for TV = 4000 and that in the box 
region in Fig. 8 for N = 64000 in the same runs. The 
four lower panels of Fig. 14 give the corresponding dis- 



placements Ari{to -f t', to + + 200) with = 10000 or 
10200, which exceed 0.3 in this short time intervals of 
200. The differences between the consecutive patterns 
are evidently larger for N = 64000 than for N = 4000. 
We can see that this system-size dependence originates 
from the large-scale vibrational motions. 

In Table I, we give the numbers of the particles with 
(a) Bi > 0, (b) Bi > and Ti = 0, (c) Bi = Ti = 0, 
and (d) Ar^ > 0.3 in the four snapshots in Fig. 14. We 
also give those of the particles commonly depicted in the 
consecutive snapshots. In Table II, the numbers of the 
B and non-B particles are given at t = 10200 for (a) 
Ar^ > 0.3 and Ti = 0, (b) Ar^ > 0.3 and = 1, and 
(c) Ar^ < 0.3 and Ti = 1. There is no particle with 
Ar^ < 0.3 and J^i = 0. We recognize the following, (i) 
About 20 — 30% of the particles with J^i = change into 
those with J^i = 1 and vice versa in a short time of 200. 
(ii) About 30% of the B particles satisfy = 1 at each 
time because of the presence of another particle j within 
their initial circles. As stated at the end of Sec.IVA, most 
of the particles with Ar^ > 0.3 and J^i = 1 are B particles 
having broken bonds (which is 97% in the example of 
Table II). (iii) About 5% of the B particles become the 
non-B particles and vice versa in a short time of 200. (iv) 
About 85% of the B particles move outside their initial 
circle to have Ar^ > 0.3. The remaining 15% B particles 
stay within their initial circle having broken bonds after 
long-distance motions of the neighboring particles. 

Dauchot et at. [5F performed an experiment on a 2D 
dense granular packing under cyclic shear near the jam- 
ming transition. Their snapshots of Ar^ and 1 — J^i {q^ 
in their notation) resemble those in Fig. 13. 



VI. THREE-DIMENSIONAL RESULTS 

Also in 3D, the four-point corelations arise from the 
bond-breakage motions and the thermal vibrational mo- 
tions. The former grow slowly with structural relax- 
ations, while the latter fluctuate relatively rapidly. 

In Fig. 15, we give snapshots of the particles at t = 
10000 and 10200 for T = 0.24 and N = 10^, where 
Ta ~ 10^. The system length is L = 23.2. Here, the 
fraction of the B particles is ^^(t) ~ 0.24 and that 
with = is 04 (t) ~ 0.15. In (a), we display the 
relatively large displacements Avi = ri{to + t) — ri{to) 
with Ari > 0.3. String-like motions are conspicuous 
[HI im [23l ISS , around which collective motions with 
lAr^l > 0.3 tend to be induced. We also display the 
particles with Ti = grouping them into B particles 
(in red) and non-B particles (in yellow) in time intervals 
[to, to + 10000] in (b) and [to, to + 10200] in (c). 

In Fig. 16, we display the B particles in (a) and (a'), 
the B particles with J^i = in (b) and (b'), and the non- 
B particles with Ti = in (c) and (c') for t = 10000 in 
the left and 10200 in the right. We use the same data as 
in Fig. 15. We can see that the B particles little change, 
but the non-B particles much change in a time interval of 
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(a) An 



[fo,fo+ioooo] (b) Ti >0 



[^o,^o+10ooo] (c) Ti>0 



[fo,fo+10200] 






FIG. 15: Snapshots in 3D for T = 0.24 and N = 10000. (a) Arrows indicate displacements Ar^ (to, to + 10000) with | Ar^l > 0.3, 
whose number is 1818. Particles with Ti — ^ classified into B particles (red) and nou-B particles (yellow) for time intervals 
[to, to + 10000] in (b) and [to, to + 10200] in (c). Depicted non-B particles are produced by the vibration modes and are 
fluctuatiing in time, while B particles are not much changed between these two times. 



(a) Bi>Q 



[fo,fo+ioooo] (a') Bi>0 



[fo,fo+10200] 




(b) ^i=0 i5^>0 [fo,fo+10000] {b') Fi=0 Bi>0 [fo,fo+10200] 




FIG. 16: Snapshots using the data in Fig. 15 for time intervals 
[to, to + 10000] (left) and [to, to + 10200] (right). Displayed are 
particles with 23^ > in (a), those with Bi > ^ and = in 
(b), and those with = J^^ = in (c). 



200. The aggregates of the B particles have grown from 
the strings in (a) in Fig. 15. The number of the total B 
particles (in the top panels) is considerably larger than 
that of the B particles with = (in the midde panels) , 
which are 2351 and 641, respectively, at t = 10000. This 
is because the particles surrounding each string can have 
broken bonds without their long-disance motions. 



Table III gives the numbers of the particles with (a) 
Bi > 0, (b) Bi > and Ti = 0, (c) Bi = Ti = 0, 
and (d) Ar^ > 0.3 in the two snapshots in Figs. 15 and 
16. Also given in the last line are those of the particles 
commonly depicted in the consecutive snapshots. On 
the other hand. Table IV presents the numbers of the B 
and non-B particles at t = 10200 for (a) Ar^ > 0.3 and 
Ti = 0, (b) Ari > 0.3 and Ti = 1, and (c) Au < 0.3 
and J^i = 1. Conspicuous features are as follows, (i) 
About 50% particles with Ar > 0.3 change into those 
with Ari < 0-3 and vice versa in a short time interval 
of 200. Only 30% of the non-B particle are common in 
the two consecutive snapshots, (ii) About 60% of the 
B particles satisfy Ar^ < 0.3 and J^i = 0. (The cor- 
responding percentage is about 15% in 2D in Table II.) 
This is because of the larger coordination number or the 
larger bonded particles around each particle in 3D. That 
is, there are about 10 bonded particles around a parti- 
cle which has participated in a string-like motion and 
moved over a molecular distance, (iii) The particles with 
Ari > 0.3 and J^i = 1 amount to 366 and their fraction in 
those with Ari > 0-3 is about 0.2. This means that the 
distinct part is negligible as compared to the self part in 
Ti in Eq.(4.2), in accord with the calculation by Lacevic 
et a/. [23]. (iv) Among 366 particles with Ar^ > 0.3 and 
J^i = 1, most of them (339) are B particles. 



13 



TABLE III: Particle numbers in Figs. 15 and 16 at t = 10000 
and 10200 for N = 10^ in 3D. 



t 


B>0 


B>0,J^ = 


B = J^ = 


Ar > 0.3 


10000 


2351 


641 


811 


1818 


10200 


2393 


675 


898 


1943 


common 


2094 


400 


268 


980 



TABLE IV: Numbers of non-B and B particles for three cat- 
egories at t = 10000 for N = 10^ in 3D in Figs. 13 and 14. 





Ar > 0.3 

j- = o 


Ar > 0.3 
J=-= 1 


Ar < 0.3 
J=-= 1 


total 


8 = 


811 


27 


6811 


7649 


B>0 


641 


339 


1371 


2351 


total 


1452 


366 


8182 


10000 



VII. SUMMARY AND REMARKS 

We have examined the dynamic heterogeneity of glassy 
particle systems in the bond-breakage scheme [I3l [14] 
and in the four-point scheme [23 . The former treats the 
irreversible configuration changes, while in the latter also 
included are the reversible particle displacements due to 
the low- frequency vibrational modes. These two kinds of 
motions are both highly heterogeneous in glassy states. 

Our main results are as follows. 

(i) In Sec. Ill, we have generalized the bond breakage 
theory |l3l [14] to define the broken bond number 
Si(to,ti), the fractions of the particles with k broken 
bonds (j)i){t^k)^ the correlation function Gbir^t)^ the 
structure factor Si){q^t)^ and the susceptibility Xb{t) 
in Eqs.(3.12), (3.14), and (3.19)-(3.21). We have de- 
fined the bond-breakage time r5 in Eq.(3.4) and the 
bond-preserving time r^p in Eq.(3.16) in addition to the 
relaxation time from Fs{q,t) in Eq.(3.6). In Fig. 3, 
Xb{t) exhibits a maximum as a function of t, yielding 
the maximization time t^^^^ while the Ornstein-Zernike 
fitting of Sbiq.t^'''') yields ^5 = Fig.4. These 
quantities are nearly independent of the system size as 
long as 1 <C ^6 <^ ^• 

(ii) In Sec. IV, we have discussed the four-point theory, 
where the overlap function w{r) in Eq.(4.3) defines 
the initial circles (spheres) in 2D (3D). The overlap 
number J-i{to^ti) in Eq.(4.2) determines the correlation 
function G4{r^t)^ the structure factor S4{q^t)^ and the 
susceptibility X4(^) Eqs.(4.12)-(4.14). Maximization 
of X4(^) with respect to t yields the characteristic time 
^max have shown that the nonoverlap motions from 
the initial circles stem from the thermal excitation of the 
low-frequency vibrational modes and the escape jumps 
from temporary cages as in Figs. 5, 6, 11, and 13. The 
thermal collective motions appear from the initial stage 
{t > 10) as in Fig. 7, while the jump motions emerge 
very slowly. The maximization procedure of X4(^) 
with respect to the overlap length [23] is to maximize 
the contribution from the thermal collective motions 
to X4(t). In 2D, the system-size dependence of the 



four-point correlations is strong at long wavelengths 
even for ^4 <C 

(iii) In Sec.V, we have compared the relaxation times r5, 
^pbj Taj tb^^^j and ^4^^^ in Fig. 12 to obtain the sequence 
(5.1), where r^i^^ t^^^) is considerably shorter than r^. 
Next we have presented snapshots of the displacements, 
Bi, and Ti at t = 10^ for N = 4000 with marked 
large-scale heterogeneities in Fig. 13. We have grouped 
the particles with J^i = into B and non-B particles, 
where they are those with and without broken bonds, 
respectively. The patterns of the B particles with = 
closely resemble those of the total B particles. The 
non-B particles with J^i = arise from the low-frequency 
vibration modes undergoing relatively rapid temporal 
variations, as can be seen in the inset of Fig. 6 and in 
Fig.l4. 

(iv) Also in 3D, the four-point correlations arise from 
the thermal collective motions with Bi = and the 
bond-breakage motions with > 1 as in Figs. 15 and 
16. As a charcteristic feature in 3D, Table IV shows 
that about 60% of the B paricles satisfy Ar^ < 0.3 and 
J^i = 1. These particles suround the particles which 
have undergone string-like motions. 

We make some remarks in the following. 

(1) The heterogeneity exhibited by the low- frequency 
vibration modes still remains largely unexplored [35H4T] . 
In future work, we should examine how it depends on 
the size ratio and the composition [45H47j. 

(2) The vibration modes determine the plateaus of the 
time-correlation function Fs{q^t) in Eq.(3.5) and the 
mean square displacement M{t) in Eq.(A4), resulting 
in significant system-size effects [56l [57] . They also give 
rise to the system-size dependence of the four-point 
structure 6*4 (g', t) at small q in (b) of Fig. 10. Our present 
analysis is mostly for 2D, but Eqs.(A2) and (A8) provide 
one possible sourse of the finite size effect in 3D. 

(3) We also comment on the effect of a thermostat, 
which was used only in preparing the initial states. We 
have found that a thermostat can strongly affect the 
low- frequency vibration modes (not shown in this paper). 
For example, the second peak of X4(^) at t = L/2c± 
for T < 0.80 in Fig. 9 disappeared in the presence of a 
thermostat, presumably because it effectively increases 
the acoustic damping. 
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Appendix: Thermal positional fluctuations and 
finite size effect 

In solids, the vibration modes give rise to the ther- 
mal fluctuations of the particle displacements Ui{t) = 
ri{t) — fi with fi being the time-averaged positions. In 
our theory, we may define fi on timescales shorter than 
the bond-preserving time Tbp in Eq.(3.16). The contribu- 
tions from the large-scale modes may be calculated using 
the classical linear elasticity theory, which should be valid 
at sufficiently long wavelengths even in glass [35^. 

In 2D, we may expresss the displacement variance 
= Ei(l^iP)/^ as in Eqs.(4.6) and (4.7). For 2D 
solids ^54], use has been made of the relation, 

{\u,{t) - Uj{t)\^) = 2Ci ln(r,,/ao), (Al) 

where Ci is given in E.(4.7), r^j is the distance be- 
tween particles i and j, and ao is a microscopic length. 
Thus, Eq.(Al) represents the anomalous long-range cor- 
relation in 2D solids. These expressions folllow if the 
discrete sums over the long-wavelength vibration modes 
i^q are replaced by the wave- number integral 

(/ dqq"^-^). In 3D, there is no long- wavelength diver- 
gence, but the lower bound of the wave number (ex L~^) 
yields the following L dependence, 

{\u\^)=Do-D^/L, (A2) 

where I^o and Di are functions of T. If we perform the 
corresponding discrete summation over the modes under 
the periodic boundary condition, we obtain 

A=0.21.,Tg + ^^^). (A3) 

For t <C Tfop, we may set ri{to) — ri{to -\-t)= Ui{to) — 
UiitQ + t) by neglecting the configuration changes. The 
mean square displacement M{t) is written as 

M{t) = Y,{\r,(t,)-ri{t,+t)\^)/N 



^ Y^{\u,{t^)-Ui{U+t)\'')/N. (A4) 

i 

If the cross correlation Xli('^i(^o + ^) • Ui{tQ))/N decays 
to zero due to the acoustic damping before the a relax- 
ation, a well-defined pleteau Mp appears in M(t) with 

M(t)-Mp = 2(|np), (A5) 

for 1 ^ t ^ as already given in Eq.(4.6). 

The time correlation function Fs{q,t) in Eq.(3.5) also 
assumes a well-defined plateau value /p = /p(T, A") at 
low T. For the Gaussian distribution of n^, we have [59] 

/p - exp[-q'Mj2d] = exp[-q^{\u\^) /d], (A6) 
If the structural relaxation time is defined by Eq.(3.6), 
its dependence on A(oc A^/'^) can arise from Eq.(4.6) in 
2D and from Eq.(A2) in 3D. 
For 2D, /p depends on A as 

/p oc oc A^-«'c?i/4. (A7) 

For T = 0.56 and g' = 27r, our numerical analysis gives 
/p = 0.6 for A = 4000 in Fig.l and /p = 0.5 for A = 
64000. The ratio of these two values 0.6/0. 5 = 1.2 is close 
to the theoretical ratio 16^^^^/^ = 1.17 from Eq.(A7). 
For 3D, Eqs.(A2), (A5), and (A6) yield 

/p(T, A)//p(T,oo) = exp[BfN-'/% (A8) 
Under the periodic boundary condition, Eq.(A3) gives 

Bf ^ O.Uq^kBTn^^^/fi, (A9) 

where K is assumed to be considerably larger than /i. 
Kim and Yamamoto [57 studied the finite size effect us- 
ine; the soft-core potential for A = 108,10^, and 10"^ 
in 3D. Their data of Fs{q,t) at t = 10 for q = 2it 
and T = 0.267 may be approximately fitted to the 
form oc exp[0.25A-^/^], while Eqs.(A8) and (A9) give 
/p oc exp[0.29A~^/^] (where we obtain /i = 5 from the 
stress-strain relation). 
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